Properties of Liquid Iron along the Melting Line up to the Earth-core Pressures 
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We report a molecular dynamics study of transport coefficients and infinite frequency shear mod- 
ulus of liquid iron at high temperatures and high pressures. We observe a simultaneous rise of 
both shear viscosity and diffusion coefficient along the melting line and estimate if liquid iron can 
vitrify under Earth-core conditions. We show that in frames of the model studied in our work iron 
demonstrates a moderate increase of viscosity along the melting line. It is also demonstrated that in 
the limit of high temperatures and high pressures the liquid iron behaves similar to the soft spheres 
system with exponent n « 4.6. 

PACS numbers: 61.20.Gy, 61.20.Ne, 64.60. Kw 
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I. INTRODUCTION 

The behavior of iron at the Earth core pressure- 
temperature conditions is a topic of hot debates for many 
decades. Apparently the behavior of iron is of great im- 
portance for understanding of the phenomena occurring 
in the inner and outer core. However, it appears to be 
very difficult to find any unambiguous information about 
iron at such extreme conditions. The problem is clear: it 
is impossible to carry out direct experiments at so high 
temperatures. Therefore one has to use extrapolation 
of the lower temperatures results. This leads to a great 
variation of the most principle data. Even the location of 
the melting line of iron at high pressures is not clear: the 
results from diamond anvil experiments give the melting 
temperature approximately twice smaller then the shock 
wave experiments (a possible explanation of this discrep- 
ancy was proposed in Ref. [l[). 

The situation with transport properties of iron at high 
pressure is even much worse. The difference between the 
viscosities obtained by different methods achieves 10 14 
. Secco classified the iron viscosity estimations in three 
groups H: the ones from geodesic measurements and 
seismological investigations give the Earth-core viscosity 
up to 10 11 Pa-s. The viscosities obtained from the Earth 
magnetic field are of the order of 2.7 • 10 7 Pa ■ s. Finally 
the theoretical predictionsgive the iron viscosities from 
2.5 • 10~ 3 up to 50 Pa ■ s [3]. Obviously the discrepancy 
of 10 14 between different methods can not be recognized 
as acceptable. 

One of the possible sources of errors in the high pres- 
sure iron viscosity estimation can originate from the ex- 
trapolation of low pressure data far beyond the range of 
pressures where we have experimental data. In this re- 
spect it is reasonable to find a system which we can study 
with reasonable precision in wide range of pressures and 
temperatures. The most obvious way to implement this 



idea is to employ some model of iron in molecular dy- 
namics. 

Several authors carried out MD simulations of iron at 
Earth-core conditions with different empirical potentials 
or by means of DFT method. In Refs. 0, Q the authors 
use the same parametrization of embedded atom poten- 
tial (EAM) for iron proposed by Sutton and Chen [||. 
The authors of Ref. [j] and Ref. [f| have chosen different 
density-temperature points which makes more difficult to 
compare theirs results. However, the viscosity data from 
both articles look to be consistent with each other. The 
viscosities at the Earth core conditions obtained in both 
cited papers are of the order of magnitude of 0.01 Pa ■ s. 
In Ref. [B| the viscosities are also compared to the ab- 
initio calculations of Alfe et. al (Table I of Ref. 0). 
One can see from this comparison that the data from clas- 
sical MD calculations and from ab-initio MD are close to 
each other and that there is no systematic deviation of 
EAM data from the ab-initio ones. 

A lot of information about iron at high pressure was 
recently obtained by ab-initio simulations. It includes 
the melting line calculations @, Q , transport properties 
of iron Q and mixtures of iron with other elements [Io| . 
However, ab-initio simulations are very computationally 
expensive and do not allow to study the properties of iron 
in a vast region of pressure - temperature points. 

The goal of the present study is to perform a system- 
atic study of liquid iron properties along the melting line 
in a wide range of temperatures and pressures. This will 
allow us to see the changes which take place in a model of 
realistic liquid under huge change in pressure and temper- 
ature. Since we use the same model and all data points 
are obtained directly we do not use any extrapolation 
procedures which allows us to see the general trends in 
the liquid iron behavior along the melting line up to the 
very high pressures. 
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II. SYSTEMS AND METHODS 

As it was stated above the goal of this article is to 
study the properties of liquid iron along the melting line. 
In the literature there is a lot of different data on the 
melting line of iron obtained by different groups and us- 
ing different methods. Several authors reported melting 
line of iron from molecular dynamic simulation. 

The most extensive simulation of melting line of iron 
was performed by Belonosho et. al [111 ] . Basing on the 
EAM potential for iron introduced in this work the au- 
thors computed the melting line from P — 60 Gpa up 
to 400 GPa. This corresponds to the temperatures up 
to 8000 K which is even above the estimated Earth core 
temperature. 



Basing on these well documented data from Ref. [11 1 
we compute the properties of iron along the melting line. 
For doing this we simulate a system of 3456 iron atoms 
in the (P,T) points located along the melting line from 
T low = 2500 K up to T hlgh = 8000 K with the tempera- 
ture step dT = 100 K. Firstly we simulate the system in 
NPT ensemble to find the equilibrium density of liquid 
at the given temperature and pressure. At this stage we 
carry out 10 6 MD steps with the step size dt = 0.001 ps. 
When the equilibrium density is found we simulate the 
system in NVT ensemble at this density in order to cal- 
culate the equilibrium structure and infinite frequency 
shear modulus of the liquid. At this stage the system 
is propagated 10 6 time steps with dt — 0.0002 ps. Fi- 
nally using the obtained final structure at the chosen den- 
sity and temperature as initial conditions we carry out 
NVE simulation of the sample to find the diffusion coef- 
ficient. For doing this we simulate the system for 2.6 ■ 10 6 
steps with dt = 0.0002 ps. The infinite frequency shear 
modulus is evaluated as Gi n f — j3V < P xy >, where 
P = l/(fcgT) and P xy is off-diagonal pressure compo- 
nent. The diffusion coefficient is computed from the slope 
of mean square displacement of the particles via Einstein 
relation. 

One of the central quantities of our analysis is the shear 
viscosity of liquid iron along the melting line. In order 
to compute the viscosity we employ the Reverse Non- 
equilibrium MD method also known as Miiller-Plathc 
method [l2| . In this method an artificial momentum flux 
is imposed in the system which results in a linear veloc- 
ity profile of particles. The viscosity coefficient can be 
calculated from the slope of the velocity profile and the 
momentum transferred to the system 12j . For the viscos- 
ity calculation the system was simulated for 2 • 10 6 time 
steps with dt = 0.0002 ps. The momentum transfer was 
undertaken every 10 steps. The temperature was held 
by coupling to Berendsen thermostat with time constant 
tB = 10 • dt. The first half of the simulation was used 
for equilibration while during the second half the veloc- 
ity distribution was written in the file every 100 steps. 



All of these distributions written to the file were used to 
estimate the viscosity. 

All simulations reported in this work were done by 
LAMMPS simulation package [Hj]. 



III. RESULTS AND DISCUSSION 

As it was mentioned in the introduction we are inter- 
esting in the behavior of the transport coefficients of liq- 
uid iron along the melting line. In our previous publica- 
tions we analyzed the behavior of simple liquids along the 
melting curve 14, ljj . Two simple models were studied 
- soft spheres (<l>(r) = e(-) n with n = 12) and Lennard- 
Jones ($(r) = e ■ ((^) 12 - (-) e )) liquids. It was shown 



that in these models both shear viscosity and diffusion 
coefficient grow up upon increasing pressure along the 
melting line. However, even if these model liquids can 
be used as a rough model to analyze the most general 
trends in liquids, the simplicity of these models can re- 
sult in large qualitative errors at high pressures and high 
temperatures comparing to the experimental liquids. In 
this respect the present work is aimed to understanding 
the high temperature - high pressure behavior of an ex- 
ample of real liquid. We choose a particular case of liquid 
iron due to its importance in geophysical investigations. 

Fig.[U(a) shows the melting line of iron [ll| and Fig.Q] 
(b) demonstrates the liquid branch of the melting line in 
the density - temperature plane. This line was obtained 
by perfor ming NPT simulations at the PT data points 
from Ref. [llj . As it is expected the density of the liquid 
quickly rises along the melting curve. 

Figs. [5] (a) and (b) represent the diffusion coefficient 
and shear viscosity of liquid iron along the melting line. 
As in the case of simple models both the diffusion coeffi- 
cient and the viscosity rapidly increase with increasing 
the temperature. The viscosity rise is especially dra- 
matic: at the highest temperature studied (8000 K) it 
is 2.5 times higher then at the lowest one (2500 K). At 
the same time the diffusion coefficient increases just 1.5 
times. Anyway the liquid becomes more viscous and 
more diffusive at the same time. 

The question of simultaneous growing up of diffusion 
coefficient and shear viscosity was raised up in our pre- 
vious publication 15j. This question is interesting in 
conjunction with glass transition. The most common cri- 
teria of glass transition states that the liquid experiences 
glass transition when its viscosity reaches some very high 
value. A typical convention is that the viscosity of glass 
transition is 10 13 Poise. However, it is implicitly assumed 
that the liquid looses its diffusivity under this viscosity 
grows. However, in case of high temperatures and high 
pressures both viscosity and diffusivity increase, so one 
comes to a contradiction to the usual common view on 
glass transition. 

In order to solve this contradiction we proposed in Ref. 
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FIG. 1: (Color online) (a) Iron melting line as obtained in Ref. 
[IH in P — T oordinates. (b) Liquid branch of the melting 
line. Squares - MD data, continuous line - soft spheres-like 
approximation (see the text); 



FIG. 2: (Color online) Diffusion (a) and shear viscosity (b) 
along the melting line. Squares - MD data, continuous line - 
soft spheres approximation (see the text). 



[la] to use one more criterium of glass transition: the liq- 
uid undergoes the glass transition if the relaxation time 
becomes as long as the typical time of experiment. Dif- 
ferent publications propose to use or 100 seconds or 1000 
seconds as the glass transition relaxation time. Following 
this definition one needs to see the behavior of the relax- 
ation time in order to understand if the liquid vitrifies 
following the melting line up to extremely high temper- 
atures - high pressures limit. 

The relaxation time can be computed via Maxwell re- 
lation 16 1 



Ginf 



(1) 



where rj is the viscosity of the liquid and Gi n f the infinite 
frequency shear modulus. In the majority of experimen- 
tal situations it is supposed that the viscosity changes 
much faster then Gi n f and the relaxation time is mainly 
determined by the viscosity behavior. This case the vis- 
cosity criterium of glass transition becomes equivalent to 
the relaxation time one. However, as it was shown in 
Ref. [l5| the infinite frequency shear modulus of liquid 
can dramatically grow along the melting line. 
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soft spheres model is especially simple since it demon- 
strates a set of scaling properties along the melting line 
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FIG. 3: (Color online) Infinite frequency shear modulus (a) 
and relaxation time (b) along the melting line. Squares - MD 
data, continuous line - soft spheres approximation (see the 
text). 



Fig. [3] (a) shows the infinite frequency shear modulus of 
iron along the melting line. One can see that Gi n f dras- 
tically increases with increasing the temperature. The 
ratio of Gi n f at the highest and the lowest temperatures 
is approximately 10, while for the viscosity it is 2.5. As a 
result in spite of the rise of viscosity the relaxation time 
still decreases (Fig. [3] (b)). 

From the results presented above one can see that the 
qualitative behavior of liquid iron at high temperatures - 
high pressures is equivalent to the behavior of simple liq- 
uids such as soft spheres and Lennard- Jones ones. The 



17H20J ■ The scaling relations for all quantities presented 



here are given in Refs 
the sake of completeness 



14 , HI . Here we repeat them for 



G f ~ 7^+3 A 



(2) 



(3) 



(4) 



(5) 



In order to see the relations between the simplest 
model studied and the current liquid iron system we fit 
all the quantities above (D, rj, P and Gin/) to the rela- 
tions of the form X = a ■ T a + b, X is the quantity of 
interest and a is the correspondent exponent from eqs. 
(2) - (5). In order to get all of the exponents consistent 
with the case of soft spheres all quantities were fitted 
simultaneously. The results of such fitting are given in 
Figs. []]- El The exponent coefficient n is found to be 
equal to n — 4.568. One can see that except the melting 
pressure all of the quantities are well represented by the 
soft spheres-like scaling relations. In the case of pressure 
the deviation does not exceed 12% in the whole range 
of temperatures considered in this work. However, the 
slope from the scaling formula and from MD data are 
very different which means that the melting line itself is 
poorly represented by the scaling law. At the same time 
the transport coefficients and elastic properties (G;„/) 
are well described by the soft spheres-like model. 

It is well known that the structure of liquid metals 
can be well approximated by simple hard spheres model 
21 1 . In Ref. [22J experimental measurements of liquid 
iron structures factors were reported and the comparison 
of experimental curves with the hard spheres model was 
done. As it follows from Fig. 3 of Ref. [22| the structure 
factors of iron can be sufficiently well represented by hard 
spheres ones which proves that the main contribution 
into the liquid structure comes from repulsive part of the 
interaction. It is well known that the structure of liquid is 
closely related to its transport properties and therefore 
one can expect that simple purely repulsive models of 
liquid can reproduce the diffusion and shear viscosity of 
iron sufficiently well. 

At the same time melting line is strongly affected by 
the presence of attractive terms in the interparticle inter- 
action potential [23| which means that a purely repulsive 
model such as soft spheres should fail to reproduce the 
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melting curve of a system with both repulsive and attrac- 
tive interactions. 

One of the experimental studies of the liquid iron vis- 
cosity at high pressure was reported in Ref . |24( . The au- 
thors of this paper measured the viscosity at temperature 
as high as 2050 K and found that the change in viscosity 
comparing to the room temperature is small. Basing on 
this result they concluded that the statement proposed 
by Poirier that the viscosity of liquid is nearly constant 
along the melting line [lH is correct. From our results 
we can conclude that it is just partially true. The viscos- 
ity change along the melting line is not fast: it increases 
2.5 times on 3.2-fold temperature change. However, we 
observe the systematic rise of viscosity along the melt- 
ing line so one can not claim that it is constant: if on 
measures the viscosities along the melting line for large 
enough temperature interval one will clearly see the rise 
of viscosity. However, the temperatures studied in our 
work range from 2500 K up to 8000 K which exceeds the 
range of temperatures reported in most of experimental 
works. It means that in the range of temperatures ex- 
plored in experiments the viscosity change can be small 
enough to use Poirier statement with sufficient accuracy. 



CONCLUSIONS 

The present article represents a molecular dynamics 
study of transport coefficients and glass transition of liq- 
uid iron in the limit of high temperatures - high pressures 
along the melting line. We show that both shear viscos- 
ity and diffusion coefficient increase along the melting 
line. However, due to very rapid increase of the infi- 
nite frequency shear modulus the relaxation time drops 
quickly with increasing the pressure along the melting 
curve which means that liquid iron becomes harder to 
vitrify at higher temperatures and higher pressures. 

It is worth to note that the magnitude of viscosities we 
obtain are consistent with other simulations of liquid iron 
at high temperatures and high pressures @, 111 0] • How- 
ever, it looks that all simulations of iron at such extreme 
conditions strongly underestimate the shear viscosity. 

Surprisingly, the behavior of liquid iron at Earth-core 
like temperatures and pressures can be sufficiently well 
qualitatively described by soft spheres model which is one 
of the simplest models of liquid. By fitting the MD data 
to the soft spheres scaling relation we find that liquid 
iron is qualitatively similar to the soft spheres with the 
exponent n — 4.568. 

The most important difference between the soft 
spheres and liquid iron represented by the EAM poten- 
11| is in the collective nature of the later. It 



tions confirm this speculation and propose that the exact 
results for soft spheres reported in our previous work 15 1 
can be extrapolated to the high temperature - high pres- 
sure limit of liquids in general. 
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tial Ref. 

is well known that in the limit of high pressures the par- 
ticles come very close to each other and the system is 
dominated by the repulsive excluded volume effects and 
the collective effects can become negligible. Our Simula- 
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